********************************************************************************
* Graph all figures 
********************************************************************************

cap restore
clear all


global output ".../Figures"
global hrsdata "...data/HRS MiCDA/Tables"
global discl "...data/FSRDC/Tables"
global othdat ".../data"
global modeloutput "...data/model-output"

*-------------------------------------------------------------------------------
*Comparison of plans
*-------------------------------------------------------------------------------

clear
set obs 21
gen mw = 0 if _n==1
replace mw = mw[_n-1] + .005 if mw==.

*TSP
gen mw_tsp = .
replace mw_tsp = .01 + mw if mw<=.03
replace mw_tsp = .04 + .5*(mw-.03) if mw>.03 & mw<.06
replace mw_tsp = .05 if mw>.05

*Avg HRS
gen mw_hrs = .
replace mw_hrs = .011 + .39*mw if mw <=.065
replace mw_hrs = .0364 if mw>.065

*Most generous common private sector plan 
gen mw_mode = .
replace mw_mode = mw if mw<=.04
replace mw_mode = .04 if mw>.04

twoway line mw_tsp mw, lw(medthick) lc(black) || ///
       line mw_hrs mw, lw(medthick) lc(red) lpattern(dash) ///
	   xtitle("Worker contribution rate") ytitle("Firm contribution rate") yscale(ra(0 .05)) ylab(0(.01).05) ///
	  legend(on order(1 "TSP" 2 "Average from HRS data")  position(5) ring(0) size(medsmall) cols(1)) graphr(c(gs16))	   
	   graph export "$output/dc-plan-comparison.pdf", replace
	   
*-------------------------------------------------------------------------------
*Time series contributions and funding ratios
*-------------------------------------------------------------------------------

import excel using "$othdat/private-pension-plan-bulletin-historical-tables-and-graphs.xlsx", sheet("17") cellrange(A8:F63)	clear   
drop if A == .
keep A F
rename A year
rename F se_tot_contrib
replace se_tot_contrib = se_tot_contrib/1000

tempfile contrib
save `contrib'

import excel using "$othdat/2021-pension-data-tables.xlsx", clear sheet("S-44") cellrange(A3:D35)
keep A D
rename A year
rename D funding_ratio

merge 1:1 year using `contrib'
keep if _merge==3

twoway (line se_tot_contrib year, lw(medthick) lpattern(solid) lc(blue)) ///
	   (line funding_ratio year, lw(medthick) lc(red) lpattern(dash) yaxis(2)), ///
	   xtitle("Year") ytitle("Dollars (billions)") ytitle("Funding ratio", axis(2)) ///
	   yscale(ra(.75 1.8) axis(2)) ylabel(.8(.2)1.8, axis(2)) ///
	   legend(on order(1 "Contributions (left axis)" 2 "Funding ratio (right axis)")  position(11) ring(0) size(medsmall) cols(1)) graphr(c(gs16))
	   graph export "$output/contrib-funding-ratio.pdf", replace

*-------------------------------------------------------------------------------
*DB wealth plot 
*-------------------------------------------------------------------------------

import delimited using "$hrsdata/compwealth_frz_w10_22.csv", clear


replace dbwealth = dbwealth/1000

keep if inrange(age,50,70)
twoway (line dbwealth age, lw(medthick) lc(black)), ///
	   xline(55, lp(shortdash) lw(thin) lc(black))  ///
	   xline(65, lp(shortdash) lw(thin) lc(black))  /// 	   
	   xlab(50(5)70) ylab(75(25)175) yscale(ra(50 175)) xscale(ra(50 70)) xtitle("Age") ytitle("Thousands of 2010 dollars") ///
	   legend(off) graphr(c(gs16))
	   graph export "$output/db_wealth_profile.pdf", replace
  
*Create hypothetical vars
local zage 55 56 60 61 65 66

foreach a of local zage{ 
	local la = `a'-1
	replace dbwealth0_frz`a' = dbwealth0_frz`a'/1000
	replace dbwealth0_frz`a' = dbwealth if age==`la'
	gen dc_additions_frz`a' = dcbal10_frz`a' - dcbal_2010
	replace dc_additions_frz`a' = dc_additions_frz`a'/1000
	gen penwealth_frz`a' =  dbwealth0_frz`a' + dc_additions_frz`a'
	replace penwealth_frz`a' =  dbwealth if age==`la'
}

twoway (line dbwealth age if age<=55, lw(medthick) lpattern(solid) lc(black)) ///
	   (line dbwealth0_frz56 age, lw(medthick) lpattern(longdash_dot) lc(black%75)) ///
	   (line penwealth_frz56 age, lw(medthick) lpattern(dash) lc(black)) ///
	   (line dbwealth age if age>=55, lw(medthick) lpattern(shortdash) lc(black%50)), ///
	   xline(55, lp(shortdash) lw(thin) lc(black)) ///
	   xlab(50(5)70) ylab(75(25)175) yscale(ra(50 175)) xscale(ra(50 70)) xtitle("Age") ytitle("Thousands of 2010 dollars") ///
	   legend(on order(2 "Without DC" 3 "With DC" 4 "Counterfactual DB")  position(11) ring(0) size(medsmall) cols(1)) graphr(c(gs16))
	   graph export "$output/db_wealth_freeze55.pdf", replace

twoway (line dbwealth age if age<=60, lw(medthick) lpattern(solid) lc(black)) ///
	   (line dbwealth0_frz61 age, lw(medthick) lpattern(longdash_dot) lc(black%75)) ///
	   (line penwealth_frz61 age, lw(medthick) lpattern(dash) lc(black)) ///
	   (line dbwealth age if age>=60, lw(medthick) lpattern(shortdash) lc(black%50)), ///
	   xline(60, lp(shortdash) lw(thin) lc(black)) ///
	   xlab(50(5)70) ylab(75(25)175) yscale(ra(50 175)) xscale(ra(50 70)) xtitle("Age") ytitle("Thousands of 2010 dollars") ///
	   legend(on order(2 "Without DC" 3 "With DC" 4 "Counterfactual DB")  position(7) ring(0) size(medsmall) cols(1)) graphr(c(gs16))
	   graph export "$output/db_wealth_freeze60.pdf", replace	   
	   
twoway (line dbwealth age if age<=65, lw(medthick) lpattern(solid) lc(black)) ///
	   (line dbwealth0_frz66 age, lw(medthick) lpattern(longdash_dot) lc(black%75)) ///
	   (line penwealth_frz66 age, lw(medthick) lpattern(dash) lc(black)) ///
	   (line dbwealth age if age>=65, lw(medthick) lpattern(shortdash) lc(black%50)), ///
	   xline(65, lp(shortdash) lw(thin) lc(black)) ///
	   xlab(50(5)70) ylab(75(25)175) yscale(ra(50 175)) xscale(ra(50 70)) xtitle("Age") ytitle("Thousands of 2010 dollars") ///
	   legend(on order(2 "Without DC" 3 "With DC" 4 "Counterfactual DB")  position(7) ring(0) size(medsmall) cols(1)) graphr(c(gs16)) 	   
	   graph export "$output/db_wealth_freeze65.pdf", replace

*-------------------------------------------------------------------------------
* Treatment effect plots 
*-------------------------------------------------------------------------------

*****************
*Employment 
*****************

import excel using "$discl/regression-tables-2-ipw.xls", sheet("lfp") clear

do "$output/prep_excel_sheet.do"

local stub "_2539 _4049 _5055 _5054 _5664 _2570"
foreach s of local stub{
	twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.04(.02) .06) yscale(ra(-.04 .06))
	   graph export "$output/emp`s'.pdf", replace

}

local s "_6570"
replace lb`s' = -.04 if lb`s'<-.04
replace ub`s' = .06 if ub`s'>.06
twoway rcap lb`s' ub`s' k if k<12, lwidth(thin) col(black%75) || ///
	   rspike lb`s' ub`s' k if k==12, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.04(.02).06) yscale(ra(-.04 .06))
graph export "$output/emp`s'.pdf", replace

*Placebo plots 
twoway rcap lb_2539 ub_2539 k, lwidth(thin) col(gray%80) || ///
	   rcap lb_4049 ub_4049 k, lwidth(thin) col(black) || ///
	   connected b_2539 k, mcolor(gray%80) msymbol(circle_hollow) lcolor(gray%70) lpattern(dash) lwidth(thick) || /// 
	   connected b_4049 k, mcolor(black) msymbol(circle_hollow) lcolor(black) lwidth(medthick)  /// 
	   yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(on order(3 "25-39" 4 "40-49") title(`"Age group"', size(medsmall))  position(2) ring(0) size(medsmall)) ///
	   xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) xlabel(-5(1)12) ylabel(-.04(.02).06) yscale(ra(-.04 .06))
	   graph export "$output/placebo_emp.pdf", replace

*Balanced panel	   
import excel using "$discl/regression-tables-4-10y-ipw.xls", sheet("lfp") clear

do "$output/prep_excel_sheet.do"

local stub "_2539 _4049 _5055 _5054 _5664 _2570"
foreach s of local stub{
	twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.06(.02).06) yscale(ra(-.06 .06))
	   graph export "$output/emp`s'_10y.pdf", replace

}	

local s "_6570"
replace lb`s' = -.06 if lb`s'<-.06
replace ub`s' = .06 if ub`s'>.06
twoway rcap lb`s' ub`s' k if k<12, lwidth(thin) col(black%75) || ///
	   rspike lb`s' ub`s' k if k==12, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.06(.02).06) yscale(ra(-.06 .06))
graph export "$output/emp`s'_10y.pdf", replace   

***********
*Unweighted
***********

import excel using "$discl/regression-tables-1-unwtd.xls", sheet("lfp") clear

do "$output/prep_excel_sheet.do"

local stub "_5055 _5664"
foreach s of local stub{
	twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.04(.02) .06) yscale(ra(-.04 .06))
	   graph export "$output/emp`s'_unwtd.pdf", replace

}

local s "_6570"
replace lb`s' = -.04 if lb`s'<-.04
replace ub`s' = .06 if ub`s'>.06
twoway rcap lb`s' ub`s' k if k<12, lwidth(thin) col(black%75) || ///
	   rspike lb`s' ub`s' k if k==12, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.04(.02).06) yscale(ra(-.04 .06))
graph export "$output/emp`s'_unwtd.pdf", replace


*************
*Gender split
************* 

import excel using "$discl/regression-tables-2-ipw.xls", sheet("56_64_gender_lfp") clear

*clean up
drop if _n<=4
rename B male
rename C female

drop if _n>36

*remove stars and brackets
local stub "male female"
foreach s of local stub{
	replace `s' = subinstr(`s',"*","",.)
	replace `s' = subinstr(`s',"(","",.)
	replace `s' = subinstr(`s',")","",.)	
}

*destring 
local stub "male female"
foreach s of local stub{
	destring `s', force replace
}

*flag point estimates and SEs
gen vtype = 1
replace vtype = 2 if A == ""
drop A 

gen id = _n
replace id = ceil(id/2)

*reshape
reshape wide male female, i(id) j(vtype)

gen k = id-6

*create 95%CIs
local stub "male female"
foreach s of local stub{
	rename `s'1 b`s'
	gen lb`s' = b`s' - `s'2*1.96
	gen ub`s' = b`s' + `s'2*1.96
	rename `s'2 se_`s'
}

twoway rcap lbmale ubmale k, lwidth(thin) col(red%70) || ///
	   rcap lbfemale ubfemale k, lwidth(thin) col(blue%60) || ///
	   connected bmale k, mcolor(red%70) msymbol(circle_hollow) lcolor(red%70) lpattern(dash) lwidth(thick) || /// 
	   connected bfemale k, mcolor(blue%60) msymbol(circle_hollow) lcolor(blue%60) lwidth(thick)  /// 
	   yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(on order(3 "Men" 4 "Women")  position(11) ring(0) size(medsmall)) ///
	   xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) xlabel(-5(1)12) ylabel(-.04(.02).08) yscale(ra(-.04 .08))
	   graph export "$output/gender_emp_effect.pdf", replace

*****************
*Retirement 	   
*****************

import excel using "$discl/regression-tables-2-ipw.xls", sheet("R") clear

do "$output/prep_excel_sheet.do"

local stub "_2539 _4049 _5055 _5054 _5664 _2570"
foreach s of local stub{
	   twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.06(.02).04) yscale(ra(-.06 .04))
	   graph export "$output/R`s'.pdf", replace

	
}

local s "_6570"
replace lb`s' = -.06 if lb`s'<-.06
replace ub`s' = .04 if ub`s'>.04
twoway rcap lb`s' ub`s' k if k<12, lwidth(thin) col(black%75) || ///
	   rspike lb`s' ub`s' k if k==12, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.06(.02).04) yscale(ra(-.06 .04))
graph export "$output/R`s'.pdf", replace
	   
*Placebo plots 
twoway rcap lb_2539 ub_2539 k, lwidth(thin) col(gray%80) || ///
	   rcap lb_4049 ub_4049 k, lwidth(thin) col(black) || ///
	   connected b_2539 k, mcolor(gray%80) msymbol(circle_hollow) lcolor(gray%70) lpattern(dash) lwidth(thick) || /// 
	   connected b_4049 k, mcolor(black) msymbol(circle_hollow) lcolor(black) lwidth(medthick)  /// 
	   yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) ///
	   xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) xlabel(-5(1)12) ylabel(-.06(.02).04) yscale(ra(-.06 .04))
	   graph export "$output/placebo_R.pdf", replace
	   
*Balanced panel	   
import excel using "$discl/regression-tables-4-10y-ipw.xls", sheet("R") clear

do "$output/prep_excel_sheet.do"

local stub "_2539 _4049 _5055 _5054 _5664 _2570"
foreach s of local stub{
	twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.06(.02).06) yscale(ra(-.06 .06))
	   graph export "$output/R`s'_10y.pdf", replace

}	

local s "_6570"
replace lb`s' = .06 if lb`s'>.06
replace ub`s' = -.06 if ub`s'<-.06
twoway rcap lb`s' ub`s' k if k<12, lwidth(thin) col(black%75) || ///
	   rspike lb`s' ub`s' k if k==12, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.06(.02).06) yscale(ra(-.06 .06))
graph export "$output/R`s'_10y.pdf", replace  	 

*************
* Unweighted
*************

import excel using "$discl/regression-tables-1-unwtd.xls", sheet("R") clear

do "$output/prep_excel_sheet.do"

local stub "_5055 _5664"
foreach s of local stub{
	   twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.06(.02).04) yscale(ra(-.06 .04))
	   graph export "$output/R`s'_unwtd.pdf", replace

	
}

local s "_6570"
replace lb`s' = -.06 if lb`s'<-.06
replace ub`s' = .04 if ub`s'>.04
twoway rcap lb`s' ub`s' k if k<12, lwidth(thin) col(black%75) || ///
	   rspike lb`s' ub`s' k if k==12, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.06(.02).04) yscale(ra(-.06 .04))
graph export "$output/R`s'_unwtd.pdf", replace
  
	   
*****************
*Log earnings 
*****************

import excel using "$discl/regression-tables-2-ipw.xls", sheet("log_earn") clear

do "$output/prep_excel_sheet.do"

local stub "_2539 _4049 _5055 _5054 _5664 _2570"
foreach s of local stub{
	   twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Log difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) 
	   graph export "$output/log_earn`s'.pdf", replace

	
}

local s "_6570"
twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) 
graph export "$output/log_earn`s'.pdf", replace
	   
*Placebo plots 
twoway rcap lb_2539 ub_2539 k, lwidth(thin) col(gray%80) || ///
	   rcap lb_4049 ub_4049 k, lwidth(thin) col(black) || ///
	   connected b_2539 k, mcolor(gray%80) msymbol(circle_hollow) lcolor(gray%70) lpattern(dash) lwidth(thick) || /// 
	   connected b_4049 k, mcolor(black) msymbol(circle_hollow) lcolor(black) lwidth(medthick)  /// 
	   yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) ///
	   xtitle("Year relative to freeze", height(6)) ytitle("Log difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) xlabel(-5(1)12)  ylabel(-.1(.1).4) yscale(ra(-.1 .4))
	   graph export "$output/placebo_log_earn.pdf", replace
	   
*Balanced panel	   
import excel using "$discl/regression-tables-4-10y-ipw.xls", sheet("log_earn") clear

do "$output/prep_excel_sheet.do"

local stub "_2539 _4049 _5055 _5054 _5664 _2570"
foreach s of local stub{
	twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Log difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) 
	   graph export "$output/log_earn`s'_10y.pdf", replace

}	

local s "_6570"
replace lb`s' = .06 if lb`s'>.06
replace ub`s' = -.06 if ub`s'<-.06
twoway rcap lb`s' ub`s' k if k<12, lwidth(thin) col(black%75) || ///
	   rspike lb`s' ub`s' k if k==12, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Log difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) 
graph export "$output/log_earn`s'_10y.pdf", replace  

***********
*Unweighted	   
***********

import excel using "$discl/regression-tables-1-unwtd.xls", sheet("log_earn") clear

do "$output/prep_excel_sheet.do"

local stub "_5055 _5664"
foreach s of local stub{
	   twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Log difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) 
	   graph export "$output/log_earn`s'_unwtd.pdf", replace

	
}

local s "_6570"
twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) 	   
graph export "$output/log_earn`s'_unwtd.pdf", replace

*****************
*EE-rate
*****************

import excel using "$discl/regression-tables-2-ipw.xls", sheet("ee_sep") clear

do "$output/prep_excel_sheet.do"

local stub "_2539 _4049 _5055 _5054 _5664 _6570 _2570"
foreach s of local stub{
	   twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.04(.02).04) yscale(ra(-.04 .04))
	   graph export "$output/ee_sep`s'.pdf", replace
}

	   
*Placebo plots 
twoway rcap lb_2539 ub_2539 k, lwidth(thin) col(gray%80) || ///
	   rcap lb_4049 ub_4049 k, lwidth(thin) col(black) || ///
       connected b_2539 k, mcolor(gray%80) msymbol(circle_hollow) lcolor(gray%70) lpattern(dash) lwidth(thick) || /// 
	   connected b_4049 k, mcolor(black) msymbol(circle_hollow) lcolor(rblack) lwidth(medthick)  /// 
	   yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) ///
	   xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) xlabel(-5(1)12) ylabel(-.04(.02).04) yscale(ra(-.04 .04))
	   graph export "$output/placebo_ee_sep.pdf", replace

*Balanced panel	   
import excel using "$discl/regression-tables-4-10y-ipw.xls", sheet("ee_sep") clear

do "$output/prep_excel_sheet.do"

local stub "_2539 _4049 _5055 _5054 _5664 _6570 _2570"
foreach s of local stub{
	twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.05(.05).1) yscale(ra(-.05 .1))
	   graph export "$output/ee_sep`s'_10y.pdf", replace

}	

*Unweighted

import excel using "$discl/regression-tables-1-unwtd.xls", sheet("ee_sep") clear

do "$output/prep_excel_sheet.do"

local stub "_5055 _5664 _6570"
foreach s of local stub{
	   twoway rcap lb`s' ub`s' k, lwidth(thin) col(black%75) || ///
	   scatter b`s' k, mc(black) m(Oh) lw(thick) xlab(-5(1)12) || ///
	   line b`s' k, lw(medthick) lc(black) yline(0, lp(shortdash) lw(thin) lc(gs8)) xline(0, lp(shortdash) lw(thin) lc(gs8))  ///
	   legend(off) xtitle("Year relative to freeze", height(6)) ytitle("Difference", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) ylabel(-.04(.02).04) yscale(ra(-.04 .04))
	   graph export "$output/ee_sep`s'_unwtd.pdf", replace
}

*-------------------------------------------------------------------------------
*HRS vs LEHD retirements
*-------------------------------------------------------------------------------

*LEHD
import delimited using "$discl/R_rates_2004cohort.csv", clear
keep if _n>=3

forvalues v = 1/7{
	destring v`v', force replace 
}
rename v1 year
rename v2 b_a_5055 //a = admin data
rename v3 se_a_5055
rename v4 b_a_5664
rename v5 se_a_5664
rename v6 b_a_6570
rename v7 se_a_6570

*Create 95%CIs
local stub "_a_5055 _a_5664 _a_6570"
foreach s of local stub{
	gen lb`s' = b`s' - se`s'*1.96
	gen ub`s' = b`s' + se`s'*1.96
}

tempfile admindat
save `admindat'

*HRS
import delimited using "$othdat/R_rates_allHRS04.csv", clear
rename _50_55_avg b_s_5055 //s = survey data
rename _50_55_se se_s_5055 
rename _56_64_avg b_s_5664
rename _56_64_se se_s_5664
rename _65_70_avg b_s_6570
rename _65_70_se se_s_6570

*Create 95%CIs
local stub "_s_5055 _s_5664 _s_6570"
foreach s of local stub{
	gen lb`s' = b`s' - se`s'*1.96
	gen ub`s' = b`s' + se`s'*1.96
}

merge 1:1 year using `admindat'

sort year

*Graphs
*50-55 with label
local s "_5055"
twoway rcap lb_a`s' ub_a`s' year, lwidth(thin) col(blue%80) || ///
	   rcap lb_s`s' ub_s`s' year, lwidth(thin) col(red%80) || ///
       connected b_a`s' year, mcolor(blue%80) msymbol(circle_hollow) lcolor(blue%80) lpattern(dash) lwidth(thick) || /// 
	   connected b_s`s' year, mcolor(red%80) msymbol(circle_hollow) lcolor(red%80) lwidth(thick)  /// 
	   legend(on order(3 "LEHD" 4 "HRS")  position(11) ring(0) size(medsmall)) ///
	   xtitle("Year", height(6)) ytitle("Rate", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) 
graph export "$output/comparative_R_rates`s'.pdf", replace
	   
	   
local stubs "_5664 _6570"
foreach s of local stubs{
twoway rcap lb_a`s' ub_a`s' year, lwidth(thin) col(blue%80) || ///
	   rcap lb_s`s' ub_s`s' year, lwidth(thin) col(red%80) || ///
       connected b_a`s' year, mcolor(blue%80) msymbol(circle_hollow) lcolor(blue%80) lpattern(dash) lwidth(thick) || /// 
	   connected b_s`s' year, mcolor(red%80) msymbol(circle_hollow) lcolor(red%80) lwidth(thick)  /// 
	   legend(off) ///
	   xtitle("Year", height(6)) ytitle("Rate", height(6)) ///
	   graphr(c(gs16)) xsize(12) ysize(10) 
	   graph export "$output/comparative_R_rates`s'.pdf", replace
}
	   	   
*-------------------------------------------------------------------------------
* Fit plot: model versus data
*-------------------------------------------------------------------------------

*Employment rate: data
import delimited using "$discl/Controllfp_levels.csv", clear
gen k = _n-6
rename v1 emprate_data
rename v2 se_emprate_data
gen lb_emprate = emprate_data - 1.96*se_emprate_data
gen ub_emprate = emprate_data + 1.96*se_emprate_data
drop if k==-5 | k==13

tempfile emprate_data
save `emprate_data'

*Employment rate: model
import delimited using "$modeloutput/E_S.csv", clear
gen k = _n-5		
rename v1 emprate_mod  

merge 1:1 k using `emprate_data'
drop _merge

tempfile emprate
save `emprate' 

twoway  rcap lb_emprate ub_emprate k, lwidth(thin) col(gray%80) || ///
	    connected emprate_data k, mcolor(gray) msymbol(square) msize(small) lcolor(gray%80) lwidth(thick) lpattern(dash) || ///
	   connected emprate_mod k, mcolor(black) msymbol(circle_hollow) msize(small) lcolor(black) lwidth(medthick) lpattern(solid) ///
	   xtitle("Period relative to freeze") ytitle("Rate") xlab(-3(3)12) xscale(ra(-4 12)) ///
	   legend(on order(2 "Data" 3 "Model")  position(2) ring(0) size(medsmall) cols(1)) graphr(c(gs16))
	   graph export "$output/model_fit_emprate.pdf", replace
	   
*Employment rate: IES = 1.5

import delimited using "$modeloutput/E_S_eis15.csv", clear
gen k = _n-5		
rename v1 emprate_mod_eis15

tempfile emprate_eis15
save `emprate_eis15' 
   
*Employment rate: IES = .75

import delimited using "$modeloutput/E_S_eis075.csv", clear
gen k = _n-5		
rename v1 emprate_mod_eis075

tempfile emprate_eis075
save `emprate_eis075' 

*Employment rate: \gamma = \phi = 0	 

import delimited using "$modeloutput/E_S_restr.csv", clear
gen k = _n-5		
rename v1 emprate_mod_restr  

tempfile emprate_restr
save `emprate_restr' 


use `emprate', clear 
merge 1:1 k using `emprate_eis15'
drop _merge 
merge 1:1 k using `emprate_eis075'
drop _merge 
merge 1:1 k using `emprate_restr' 


twoway  rcap lb_emprate ub_emprate k, lwidth(thin) col(gray%80) || ///
	    connected emprate_data k, mcolor(gray) msymbol(square) msize(small) lcolor(gray%80) lwidth(medthick) lpattern(dash) || ///
	   connected emprate_mod k, mcolor(red) msymbol(circle) msize(small) lcolor(red%80) lwidth(medthick) lpattern(solid) || ///
	   connected emprate_mod_eis15 k, mcolor(black) msymbol(triangle) msize(small) lcolor(black) lwidth(medthin) lpattern(solid) || ///
	   connected emprate_mod_eis075 k, mcolor(blue) msymbol(triangle) msize(small) lcolor(blue%80) lwidth(medthin) lpattern(dash) || ///	
	   connected emprate_mod_restr k, mcolor(orange) msymbol(square) msize(small) lcolor(orange%80) lwidth(medthin) lpattern(dash_dot) ///		   
	   xtitle("Period relative to freeze") ytitle("Rate") xlab(-3(3)12) xscale(ra(-4 12)) ///
	   legend(on order(2 "Data" 3 "Baseline model" 4 "IES = 1.5" 5 "IES = 0.75" 6 "{&gamma}={&phi}=0")  position(2) ring(0) size(medsmall) cols(1)) graphr(c(gs16))
	   graph export "$output/model_fit_emprate_restrictions.pdf", replace


*Treatment effects: data
import delimited using "$discl/te_emp_5664.csv", clear
drop if _n==1
rename v1 b_data
rename v2 se_b_data
destring b_data, force replace
gen lb_b_data = b_data - 1.96*se_b_data
gen ub_b_data = b_data + 1.96*se_b_data
gen k = _n-6
keep if k>=0

tempfile teffects_data
save `teffects_data'

*Treatment effects: model 
import delimited using "$modeloutput/TE_S.csv", clear
gen k = _n-1
rename v1 b_model 

merge 1:1 k using `teffects_data'
drop _merge 

tempfile teffects 
save `teffects' 

twoway rcap lb_b_data ub_b_data k, lwidth(thin) col(gray%80) || ///
	   connected b_data k, mcolor(gray) msymbol(square) msize(small) lcolor(gray%80) lwidth(medthick) lpattern(dash) || ///
	   connected b_model k, mcolor(black) msymbol(circle_hollow) msize(small) lcolor(black) lwidth(medthick) lpattern(solid) ///
	   yline(0, lp(dash) lw(thin) lc(black)) ///
	   xtitle("Period relative to freeze") ytitle("Rate")  xlab(0(3)12) xscale(ra(0 12)) ///
	   legend(on order(2 "Data" 3 "Model")  position(5) ring(0) size(medsmall) cols(1)) graphr(c(gs16))
	   graph export "$output/model_fit_teffect.pdf", replace
	   
*Treatment effects: IES = 1.5

import delimited using "$modeloutput/TE_S_eis15.csv", clear
gen k = _n-1		
rename v1 b_eis15

tempfile te_eis15
save `te_eis15' 
   
*Treatment effects: IES = .75

import delimited using "$modeloutput/TE_S_eis075.csv", clear
gen k = _n-1		
rename v1 b_eis075

tempfile te_eis075
save `te_eis075' 

*Treatment effects: \gamma = \phi = 0	 

import delimited using "$modeloutput/TE_S_restr.csv", clear
gen k = _n-1		
rename v1 b_restr

tempfile te_restr
save `te_restr' 

use `teffects', clear 
merge 1:1 k using `te_eis15'
drop _merge 
merge 1:1 k using `te_eis075'
drop _merge 
merge 1:1 k using `te_restr' 

twoway rcap lb_b_data ub_b_data k, lwidth(thin) col(gray%80) || ///
	   connected b_data k, mcolor(gray) msymbol(square) msize(small) lcolor(gray%80) lwidth(medthick) lpattern(dash) || ///
	   connected b_model k, mcolor(red) msymbol(circle) msize(small) lcolor(red%80) lwidth(medthick) lpattern(solid) || ///
	   connected b_eis15 k, mcolor(black) msymbol(triangle) msize(small) lcolor(black) lwidth(medthin) lpattern(solid) || ///
	   connected b_eis075 k, mcolor(blue) msymbol(triangle) msize(small) lcolor(blue%80) lwidth(medthin) lpattern(dash) || ///	
	   connected b_restr k, mcolor(orange) msymbol(square) msize(small) lcolor(orange%80) lwidth(medthin) lpattern(dash_dot) ///		   
	   yline(0, lp(dash) lw(thin) lc(black)) ///
	   xtitle("Period relative to freeze") ytitle("Rate")  xlab(0(3)12) xscale(ra(0 12)) ///
	   legend(on order(2 "Data" 3 "Baseline model" 4 "IES = 1.5" 5 "IES = 0.75" 6 "{&gamma}={&phi}=0")  position(5) ring(0) size(medsmall) cols(1)) graphr(c(gs16))
	   graph export "$output/model_fit_teffect_restrictions.pdf", replace
	   
